function dy = master_function_v4(t, y, x, transitions_matrix, force_rate, variables_indexes,  initial_rates)

    s = 1;
    N = length(initial_rates)/2;
    rate_matrix = transitions_matrix;
    for i = 1:length(transitions_matrix)
        for j = 1:length(transitions_matrix)
            if(transitions_matrix(i, j) > 0)
                index = find(variables_indexes==s); 
                if(sum(index)>0)
                    k_0 = x(index);
                else
                    k_0 = initial_rates(s);
                end
                
                index = find(variables_indexes==(s+N)); 
                if(sum(index)>0)
                    x_0 = x(index);
                else
                    x_0 = initial_rates(s+N);
                end
                
                %x_0
                rate_matrix(i, j) = k_0(1) * exp(t*force_rate*x_0(1)/4.1);
                
                if(rate_matrix(i, j) >= 5000)
                    rate_matrix(i, j) = 5000;
                end
                s = s + 1;
            end
        end
    end
    
    diagonal = - sum(rate_matrix);
    
    for i = 1:length(transitions_matrix)
        rate_matrix(i, i) = diagonal(i);
    end
    
    dy = rate_matrix * y;
end

